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We study transport properties of a quantum dot coupled to a Majorana zero mode and two 
normal leads. We investigate the full counting statistics of charge tunneling events which allows 
one to extract complete information about current fluctuations. Using a Keldysh path-integral 
approach, we compute the cumulant generating function. We first consider a noninteracting spinless 
regime, and find that for the symmetric dot-lead couplings, the zero-frequency cumulants exhibit a 
universal pattern of Euler numbers, independent of the microscopic parameters. For a spinful case, 
the Coulomb interaction effects are discussed for both strong interaction (single-electron occupancy 
regime) and weak interactions (perturbative regime). Compared to the case without Majorana 
coupling, we show that, while the tunneling conductance might exhibit zero-bias anomaly, the full 
counting statistics is qualitatively different in the presence of the Majorana coupling. 

PACS numbers: 73.21.Hb, 71.10.Pm, 74.78.Fk, 72.70.+m 


I. INTRODUCTION 

Majorana zero-energy modes (MZMs) have recently at¬ 
tracted enormous theoretical and experimental attention 
[1-4] due to their exotic non-Abelian braiding statistics 
[5-7] and potential application to fault-tolerant topolog¬ 
ical quantum computation [8]. A large number of theo¬ 
retical proposals has been put forward to realize MZMs 
in topological superconductors (TSCs) [9-17], see also 
reviews [ , 18, 19] for more details. One of the most 
promising proposals involves a semiconductor nanowire 
with strong spin-orbit interaction coupled to a conven¬ 
tional s-wave superconductor [13, 1- ]. An appropriate 
combination of the spin-orbit coupling, Zeeman splitting 
and induced s-wave pairing allows one to realize an effec¬ 
tively spinless p -wave superconductivity at the interface 
which is characterized by the presence of MZMs bound 
to certain defects (i.e vortices in 2D and domain walls 
in ID) [7, 20]. The simplest way to detect MZMs is 
to measure local density of states at the defect and to 
probe the emergence of the zero-energy resonance across 
the topological phase transition. The first Majorana 
tunneling spectroscopy experiment, based on a semicon¬ 
ductor/superconductor heterostructure proposal [13, 14], 
was performed in Delft [21]. Later on, the observation 
of zero bias peak in a finite magnetic field, consistent 
with the theoretical predictions [22], was reported by 
many other experimental groups [21, 23-28]. The main 
challenge of these measurements is to exclude the other 
false-positive contributions to the zero-bias peak that are 
ubiquitous in condensed matter systems such as Kondo 
effect [29, 30], disorder in the topological region[31-35] 
and in the leads [36, 37] as well as some other resonant 
Andreev scattering phenomena [19]. The feature distin¬ 
guishing the Majorana origin of the zero-bias peak from 
the other mechanisms is the quantized zero-bias peak 
conductance of 2 e 2 /h which is a universal property of 
Majorana zero modes [38, 39]. However, due to the large 
subgap conductance (so-called “soft gap” problem) ob¬ 


served in tunneling experiments [21, 23-28], the largest 
observed height of the zero-bias peak was at most 30% of 
the predicted value. Therefore, additional experimental 
tests [37, 40-55] are necessary in order to conclusively 
confirm the presence of MZMs in the semiconductor- 
superconductor heterostructures. 

In this paper, we study current correlations in a meso¬ 
scopic device consisting of a quantum dot (QD) coupled 
to a MZM and two normal leads. The possibility to tune 
the couplings between QD and other conductors as well 
as QD gate voltage allows one to study current corre¬ 
lations in a well-controlled environment. We show here 
that in the case of a symmetric left-right normal metal 
coupling, see Fig.l for a layout of the proposed device, 
current fluctuations are characterized by a universal pat¬ 
tern of the zero-frequency cumulants. We argue that the 
measurement of such cumulants allows one to exclude 
other false-positive signatures in tunneling transport and 
uniquely identify the presence of the putative Majorana 
modes. 

The transport properties of a strongly interacting QD 
coupled to a MZM and a single normal lead (NL) have 
been studied using master equations, valid in the high- 
temperature regime, in Ref. . The low-temperature 
behavior of the MZM-QD-NL system and the interplay 
between Kondo and Majorana couplings was considered 
in Ref. 43 finding that zero bias tunneling conductance 
exhibits strong temperature dependence, which is dis¬ 
tinct from that of a MZM-NL structure [38, 39]. Later 
on, Cheng, et al. [46] revisited the the low-temperature 
behavior of the MZM-QD-NL system, and found that 
Majorana coupling significantly modifies the low-energy 
properties of the QD and drives the system to a new (dif¬ 
ferent from Kondo) infrared fixed point. They also con¬ 
firmed that the temperature dependence of the zero bias 
conductance at the particle-hole symmetric point is sim¬ 
ilar to that of the MZM-NL structure [39]. The zero-bias 
conductance of a noninteracting QD with a side-coupled 
MZM (ungrounded TSC) through two normal leads was 
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considered in Refs. 42 and i5, where it was predicted 
that the tunneling conductance is given by e 2 /2 h for 
symmetric QD-lead couplings. Ref. 44 considered a spin¬ 
ful QD in the Kondo regime for this two-lead structure, 
and studied the QD spectrum and zero-bias conductance 
by using numerical renormalization group method, which 
shows that the zero-bias conductance is 3e 2 /2/i for small 
QD-MZM coupling. The shot noise of a different two-lead 
structure (with a grounded TSC) has been studied in Ref. 
5G for both noninteracting spinless QD and spinful Kondo 
QD predicting that the shot noise not only shows univer¬ 
sal behaviors but also can be used for qualitatively distin¬ 
guishing MZMs with other modes. There has been also 
an experimental interest in QD-superconductor devices. 
The interplay of the Kondo effect and superconductivity 
has been revisited in Refs. [28, 29, 57, 58]. A natural 
realization for the proposed experimental setup, see Fig. 
1, involves a T-junction of the semiconductor nanowires 
which can be grown using vapour - liquid - solid growth 
technique [59]. The QD can be created near the junction, 
two normal leads and the TSC are connected to each leg 
of the T-junction. Thus, we believe that the setup we 
propose in the paper is within the experimental reach. 

Although conductance and shot noise exhibit peculiar 
universal dependence due to the MZM coupling, it is in¬ 
sightful to obtain the full probability distribution func¬ 
tion of the charge transferred through the QD, which can 
serve as the Majorana sensor. The theory of full count¬ 
ing statistics (FCS) [60-62] for charge transport in meso¬ 
scopic systems was established by analyzing nonequilib¬ 
rium transport. A great effort has been made to investi¬ 
gate various aspects of FCS in a variety of systems the¬ 
oretically [63-76] and experimentally [77-84]. Recently, 
the FCS calculation has also been considered for electron 
transport through multiterminal networks of MZMs [75]. 
In this paper, we study FCS of charge tunneling through 
a QD with a side-coupled TSC, or equivalently a QD 
coupled to a MZM. The charge transport is measured 
between two normal leads. Here we assume that TSC is 
grounded so there is also Andreev current between (left) 
lead and the superconductor. Using the Keldysh path- 
integral approach [? ], we compute the cumulant generat¬ 
ing function. We first consider a noninteracting spinless 
QD, and find that for the symmetric dot-lead couplings, 
the zero-frequency cumulants exhibit a universal pattern 
described by Euler polynomial. This result is indepen¬ 
dent of the microscopic parameters such as QD energy 
level and QD-MZM coupling. For a spinful QD with 
a small QD Zeeman splitting, we compute FCS in the 
regime of weak (perturbative regime) and strong (single¬ 
electron occupancy regime) Coulomb interactions. In the 
former case, we compute the interaction-induced correc¬ 
tion to the cumulant generating function up to the lead¬ 
ing order in U (i.e. U 2 ). In the latter case, we apply 
a slave boson mean field approach, and study the FCS 
due to the interplay between the Kondo and Majorana 
coupling. 

The paper is organized as follows. In Sec. II, we re- 



FIG. 1. Proposed experimental setup to measure the distri¬ 
bution function of the transmitted charge, namely the full 
counting statistics. The QD is created near the center of the 
semiconductor nanowire T-junction, two normal metal leads 
are attached to the upper and lower legs, and an s-wave super¬ 
conductor is in proximity to the third lead. The latter realizes 
topological superconductor hosting two MZMs 71 and 72 . An¬ 
other electron channel capacitively couples the T—junction, 
and can be used as a charge sensor to measure the charge 
distribution function. 

view the formalism of the FCS calculation for the meso¬ 
scopic transport problem. In Sec. Ill, we introduce the 
QD-MZM model and compute the cumulant generating 
function for this model with noninteracting spinless QD. 
In Sec. IV, we consider weak Coulomb interaction effect 
for a spinful QD, and compute the leading order inter¬ 
action correction to the cumulant generating function by 
using a diagrammatic perturbation method. In Sec. V, 
we consider strong Coulomb interaction effect for a spin¬ 
ful QD, and study how Kondo and Majorana couplings 
affect the FCS within a slave boson mean field approach. 
Finally, the conclusions are shown in the Sec. VI. 


II. FULL COUNTING STATISTICS: GENERAL 
FORMALISM 

In this section, we will review the formalism for calcu¬ 
lating full counting statistics (FCS) of charge fluctuations 
in a mesoscopic system, we refer a reader to Ref. [86] for 
more details. Consider the distribution function V q for q 
electrons to be transferred through a mesoscopic device 
within the measurement time T. Here we assume that 
the measurement time is long enough (T e/I) so that 
the average number of electrons M transferred within T 
is large, i.e. M 1. The distribution function V q al¬ 
lows one to extract more information about the nature of 
the charge carriers as well as the statistics of the charge 
fluctuations. In particular, tails of the distribution con¬ 
tain information about the statistics of rare events. From 
the theoretical point of view, rather than V qi it is more 
convenient to compute the cumulant generating function 
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(CGF) x(A), defined by a Fourier transform 

X(A) = 5>^7V (1) 
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Here the auxiliary variable A represents a counting field. 
From the CGF, one can calculate the cumulants ((< 5 n q)) 
(irreducible moments of V q ) 


((S n q)) 


d n 

d(i\) n 


lny(A) 


A=o' 


( 2 ) 


for this setup (assuming that the Zeeman splitting is very 
large), and relegate the discussion of the spinful case to 
the next sections. In the former case, the model is essen¬ 
tially noninteracting, and one can calculate FCS exactly. 
The QD is coupled to two normal spinless leads which 
can be used for transport measurements. Given that 
each lead couple to the QD at a single point, one can per¬ 
form unfolding transformation and reduce the problem to 
the one corresponding to a quantum impurity coupled to 
one dimensional free fermions. Thus, the corresponding 
Hamiltonian reads 


Thus, the average current and zero-frequency sym¬ 
metrized current noise can be obtained by a simple dif¬ 
ferentiation: 


r T 


I=- dt(I(t)) = 


—i cllny 


r 


r d\ 


A=0 


(•T 


S = 


dt(6i(t)Si{ 0 ) + 6I(0)6I(t)) 


(3) 


— 1 <9 2 lny 

~T d \ 2 


A=0 

(4) 


H = -ffLead + ^QD-MZM + Ht, (7) 

where the Hamiltonians for the leads, QD-MZM system, 
the Lead-QD couplings are respectively given by 

-ffLead = -»«FEa=L,Ji/ dx^ a (x)d x lj) a (x), (8) 

#qd-mzm = e d Sd + in(d + cZ f ) 71 + *£7172, (9) 

Ht = E a =L,R (ta^(O)d + £dtya(0)) (10) 


The third cumulant and the fourth cumulant describe the 
asymmetry (or skewness) and the kurtosis (or sharpness) 
of the distribution function, and can be also straightfor¬ 
wardly obtained using x(A). 

In order to define the CGF in a proper way, Levitov 
and Lesovik introduced a “gedanken scheme”, in which 
a spin- 1/2 system magnetically couples to the electric 
current [60, 6 ]. Based on their definition, the CGF in a 
Keldysh formalism is given by 

X(A) = (t c exp [ - i J dtH x (t)] ^ (5) 

with 

H x (t)=H+^-I. (6) 

Here, the integration is preformed along the Keldysh 
contour C, Tc is time ordering; (■ • •) is the quantum- 
mechanical average [ 15]. The second term in Eq.( 6 ) 
describes the interaction between electron current and 
counting field, and A(t) has different sign on the two 
branches of the Keldysh contour, i.e. X(t) = A *4) for for¬ 
ward branch and A (t) = A *- 2 - 1 for backward branch with 
Ad) = -AC 2 ) = A. 


III. FCS FOR A SPINLESS QD WITH 
SIDE-COUPLED MZM 

A. Theoretical Model 

We now consider the setup shown in Fig. 1 - a quantum 
dot (QD) is coupled to an end of a one-dimensional (ID) 
topological superconductor (TSC) hosting two Majorana 
zero modes (MZMs) 71 and 72 at the opposite ends. For 
pedagogical reasons, we first consider the spinless model 


Here 'i/d (cd) is creation operator for an electron in the 
a-lead (QD), ed is the energy level in the QD, and t a k 
(k) is the tunnel coupling between the leads (MZM) and 
the QD. The effective Hamiltonian for the TSC is given 
in terms of the low-energy degrees of freedom (MZMs) 
assuming that the induced superconducting gap A is the 
largest energy scale in the problem. For a finite-length L 
TSC, the coupling 6 between two MFs is exponentially 
small S ~ Aexp(—L/£) with the coherence length being 
£ = vf/ A. 

We now derive the CGF for our QD-MZM model as¬ 
suming the symmetric source-drain bias (hl = eV/2 and 
Hr = —eV/2). The Hamiltonian including the counting 
field can be written as 

H X (t) = H+ (11) 

oc=L,R 

where the current operator for the a-junction is I a = 
ie[H,N a ] with the electron number operator N a for the 
a lead. One can apply a gauge transformation to remove 
the last term in the H x (t ), and obtain 

H X = Hhea.d + HqD -MZM + H x , (12) 

Ht = (t a e~ lXo ‘d ) / 2 - l pl [ (0)d+ h.c.y (13) 

a—L.R 

We note that for n = 0, the gauge symmetry of the 
Hamiltonian allows one to gauge away one of the count¬ 
ing fields so it is enough to keep the counting field in one 
of the junction. In the general case (i.e k 7 ^ 0), however, 
we need to keep both counting fields A a . 

We can now compute the path integral for the effective 
action defined by the Hamiltonian (12). Given that the 
presence of superconductor (i.e. MZM coupling) breaks 
particle number conservation, the QD Green function 
contains anomalous contributions, e.g. ( Tcd(t)d(t')} 7 ^ 
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0. Therefore, we introduce Nambu spinors: = 

(V4,V’a)/v / 2 and = {d\d)/y/2 , where a = L,i?. is 

the lead index. The effective Keldysh action now reads 


S = Steads + 'S'qD-MZM + <St, 


where 


SLe, ds =J2 I [ dtdt'& a (t)QuX(t,t')$ a (t') 

Q, J C J C 

-S'qD-MZM = [ [ dtdt'$\(t) Qq j d (t, t') (15) 

S T = - 55 J dt^tae -1 ^ 1 ip^d + c.c^j (16) 

= ~Y,J dt{^ a {t)M T} J d {t) + h.c.), (17) 


are the actions for leads, QD, and Lead-QD coupling, 
and 


Mr, a = 


t a e 


■■ >a(<) 


0 


Aa(t) 


-tte 1 — 


(18) 


Here we have already integrated out the bulk degrees of 
freedom in the leads and kept only the field (f) at the 
x = 0, i.e. at the QD. The free lead Green’s function 
Qo : a at x = 0 in the Nambu space N can be written as 


fl&M 0 


Qo.aM - ( ■ “ Q £0( w ) 


(19) 


where g^(t — t') is the P-H conjugation of g^(t — t'). 
We perform Larkin-Ovchinnikov (L-O) rotation, and the 
Green function in Keldysh space becomes 


5°M = -inpF 


= - inp F 


1 2(1-2 n a ) 
0 -1 

1 2(1-2 n a ) 
0 -1 


( 20 ) 

( 21 ) 


One notices that g% R { w) = —g°’ A (—w) and g^ K (ui) = 
~ffa' K (~ UJ )- Here, n a is the Fermi distribution function 
of the a lead with chemical potential p a , and n a cor¬ 
responds to the Fermi distribution function with —fi a . 
Assuming the symmetric source-drain bias pl = eK/2 
and pr■ = — eV/2, one can relate the Fermi function for 
particles and holes Ur = Ur and Ur = Ur. The free QD 
Green function (with MZM coupling) can be written as 


Qo,dd(^') — 


(G« dd 

^0 ,dd 

tt'R 
r 0 ,dd 

F 0 K dd\ 

0 

/~iA 
^0 ,dd 

0 

F(y dd 

F^-- 

F K TT 

G r - 

G K j j 

1 0 ,dd 

1 0 ,dd 

^0 ,dd 

,dd 

V 0 

F a — 

0 ,dd 

0 

C A i 
^0 ,dd/ 


( 22 ) 


where the retarded components read [41, 56] 


= 


( 14 ) G*-» = 


0 ,dd 
R 


to + iys + e d - SmM 

(w + ir/s ~ 2 E m (w))(w + ig s ) - e; 

u + igs ~ £d - S m (m) 

(w + irjs ~ 2 E m (w))(w + ig s ) - e; 

-E m (w) 


(23) 

(24) 




ig s - 2 E m (w))(w + iys) - ^' 
(25) 

Here the the self-energy due to MZM coupling is 
SmM = k 2 w/(w 2 — <5 2 ), and the infinitesimal gs —► 0. 
The Keldysh components are proportional to gs , and, 
thus, can be set to zero. After L-0 rotation, the action 
for the tunneling part becomes 


S t = -55 


dt 


i=cl,q 

+^( 55 

i=cl,q 


where 


M$ a = 


K, a = 


-,»« , -,aL 2 > 
e Q +e Q 






(26) 

(27) 

(28) 


are written in the Nambu space whereas g cl = I and 
7 9 = <7 1 represent the Keldysh space. Note the rela¬ 
tionship Aq 1 " 1 = — Aq 2 ^ = A q which allows one to simplify 
the expressions. After some manipulations, the cumulant 
generating function can be written as 


X(A) = J D[d t ,d]D[^,^ Q ]e i(SLo ' lds+SQD - MZM+ST) . 


Next, we perform Gaussian integration to find 


(29) 




det 


det 


^4x4 ~ Qo,ddQM\ 


14x4 — Qo,ddQMO} 


where 


(30) 


Qmx = E(E^87^0,455 M T,a^ l )\- (31) 

d i i 

Eqs.(30) is a general expression for the CGF. In the next 
sections, we will explicitly evaluate x(^) f° r different lim¬ 
iting cases. 


B. Results and Discussions 

1. QD coupled to two normal leads. 

It is instructive to review first a simple case of a non¬ 
interacting QD coupled to two normal leads. Taking the 
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limit k = 0 in Eq. (22) and substituting it into Eq. (30), 
one obtains 


r-T~ f‘OC 7 

lnX= ¥ J 00 2^ ln[(1 + T+)(1 + T " )] (32) 


with the functions T± being defined as 


well-known result for the shot noise in a QD: 

5ll _ e 2 2 d 2 lnx(A L , Xr = 0) 
eV h d\ 2 L A t =o 

_ 2e 2 4r L r fl [(r L - r R ) 2 + £ 2 ] 

" h [e 2 + (r L + r fl ) 2 ] 2 ' 1 j 

One can notice that the shot noise (as well as other cumu- 
lants) generically depend on the microscopic parameters 
for the QD such as, for example, e d . Furthermore, in the 
resonant case corresponding to e d = 0 and T/, = T R , the 
first eight cumulants are given by 

{C'r(O), • • • C 8 (0)} = {1,0,0,0,0,0,0,0}, (35) 


with C'n(O) being defined as 


= 4r L r fl 

(w ± e d ) 2 + (r L + r R ) 2 
x [n L { 1 - n R )(e i(XL ~ Xlt) - 1) + R O L\. 


(33) 


The term in the second bracket of the logarithm function 
in Eq.(32) is the particle-hole conjugation of the term in 
the first bracket. Since we consider a symmetric source- 
drain bias (eVt = —eV R = eV/2), the transformation 
w —► — co (e.g. for the terms in the second bracket) will 
result in the following changes: n R — > 1 — n R and n R —> 
1 — riL. Thus, one can see that Eq.(32) is consistent with 
the results of Ref. 68. Indeed, at zero temperature and 
to the linear order in applied bias eV, one obtains the 


C„(0) 


«*"g)) 

M 


(36) 


Notice that at the symmetric point T# = T^, all higher 
order (n > 1) cumulants become zero as e d —t 0. This 
dependence on e d is a generic feature because density 
of states in QD strongly depends on the gate voltage 
controlling e d . As we show below, this is not the case 
when QD is coupled to a TSC. 


2. QD coupled to two normal leads and a TSC 

Let us now consider a QD coupled to a TSC through 
MZM coupling, i.e. k / 0. Substituting Eqs. (19), 
(22), (27), and (28) into Eq. (30), we find the following 
expression for the cumulant generating function: 


In x(A) = 


T 


du> 
27r 


In 


1- 


Ci 


c 2 


K(A = o) ni(1 “ ni) - ic(X : bo)"' !(1 - " fi) + 


K(A = 0) 


n R {1 - n L ) + 


F 


K(A = 0) 


n L n _ R (1 - n L )(l - n R ) + 


K(A = 0) 

1 


K(A = 0) 


n L {1 - n R ) 

n L n R {n L - n R ) . (37) 


The coefficients (Ci, C 2 , Bi, B 2 , F, 1, K) in Eq.(37) 
are defined in the Appendix-A. Above expression can be 
simplified in the zero temperature limit where the terms 
proportional to n L (l-n L ), n R {l — n R ), n L n R (l-n L )(l- 
n R ), and n,Ln R (nL — n R ) vanish, and the corresponding 
expression for ln%(A) becomes 


tir)) 1 (n R (l-n L )Y = 0 (if i,j ± 0). Assuming fi L > fi R , 
the generating function can be further simplified to 


lnx(A) 


r 


did 


t-ao 'A J_ x 2tt 


In 


1 + Nitil(1 - n R ) 


N 2 n fl (l - n L ) 


(38) 


where Ni(w) = IBb/K^A = 0) and N 2 (w) = B 2 /K(A = 
0). In addition, at T = 0 we have {ul{ 1 — n R))‘ = 
n L ( 1 - n R ), (n R ( 1 - n L ))‘ = n R { 1 - n L ), and (n L { 1 - 
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lnx(A) 


T-v 0 


r 

2 


dui 
27r 


In 



A^M \ 

p( w ); • 


(39) 


where the functions are A/"(w) and T>(lo ) are defined as 
A7(w) = e -2iA H |8 e i(A i+ A R ) ri r fl [ e 2 + (EmM - W) 2 
+ (r L - r fl ) 2 ] + 4 (s M (w) 2 + 4e 2iAt r|) 

(40) 

- 4e 2lA * [ (l —e 2iAi ) E M (w) 2 r| + 2r|r fi 

+^M(w) 2 r^ + 2r L r i? (e^ + (Em(u) — w) 2 + r|j) ] |, 

T(u>) = \e 2 d + (2Em(w) — ui)u ;] + (r^, + T^) 2 (41) 

+ (r L + r H ) [2 ( e d + 2Em(w) 2 — 2Em(w)w + w“) ]. 


right-left leads is dlR/dV = e 2 / 2 h which is consistent 
with the previous work on QD coupled to an ungrounded 
TSC [42, 45]. In addition, if we reverse the right lead 
voltage Vr = —V/2 —> 14/2, the linear conductance is 
also equal to the universal value 2 e 2 //i which is expected 
based on the RG analysis [38, 4 ]. 

We now discuss higher order cumulants n > 1. One can 
check that the expressions for the shot noise (as well as 
other higher order cumulants) through the left and right 
leads are the same in the eU —► 0 limit. Beyond eV —» 0 
limit, this relation only holds for the symmetric couplings 
= Tr. Therefore, we set A r = 0 from now on and 
study current fluctuations through the left junction only. 
By expanding the CGF in terms of e lXh , one finds 




n —0 


3©’ 


AnXr 


(46) 


Eq.(39) is the main result of this section which allows one 
to compute cumulants as a function of various physical 
parameters. We now simplify above expression in the 
limit 5 = 0 and small bias eV —> 0. We keep only the 
leading order terms in eV, and simply set ui = 0 in the 
integrand. [87]. After the simplifications, we arrive at a 
very simple expression for the CGF: 


The probability V q can be obtained by the Fourier trans¬ 
form 


i r 27r 

V q = --J d\ L e~^x^L) 

" \ q)\T R + T L ) Vf^J ‘ 


(47) 


lny 

M 


= In 


eV-K) 


(r L e lXL +T R e~ iX «\ 

V r L + r H )' 


As expected, the generating function in the presence of 
(42) MZM coupling k ^ 0 is still described by the binomial 
distribution. However, the cumulants, defined as 


where M = TV/2-k = TVe 2 /h is the number of incoming 
particles during the waiting time. As mentioned above, 
the expression (71) only depends on the ratio of Tr/Tr, 
and is independent of many other microscopic parameters 
such as e d and n. This universality is due to the finite 
density of states at zero energy induced by the Majorana 
leaking into the QD, and is a characteristic feature of 
topological superconductivity. 

To get some insight we compute now currents through 
the left and right junctions. Using Eq.(3), one finds 


II 

Ir 


--- —V 

h Tr + Tr 

e 2 T r 

h Tr+Tr 


V 


(43) 

(44) 


Clearly, when F^ 7 ^ F r, there is Andreev contribution to 
the current due to the presence of a grounded supercon¬ 
ductor: 


Ia 


e 2 T r - T L 

h Tr + Tr 


V 


(45) 


C„( 0 ) 


«*"g)) 

M 


H) r 


1 d n 
MdV[ 


lnx(A) 


Al,r=0 


(48) 


follow a peculiar pattern at Tr = T# 

{C'r(O) • • • C'g(O)} = {^.0,-|,0,1,0, —} , 

Cn( 0) = (49) 

with E n {x) being the Euler polynomial. Contrary to the 
case without MZM, higher order cumulants are nonzero 
at Fi = r# and are independent of ed and k. 

The dependence of the cumulants on Majorana split¬ 
ting energy S and finite voltage bias (beyond linear in V 
contributions) can be obtained using Eq. (39). In this 
case, one needs to perform an integration of the cumulant 
spectrum C n (uj) over u € (— eV/2, eV/2), where 


C n (0j) 


((*"g))M 

M 


HY 


,1 jr 

2d\ n r 


In 1 + 


A7(w) 

DM 


Al,h=0 

(50) 


One can notice that when r# = 0, we recover the pre¬ 
vious results [46] corresponding to a single lead coupled 
to a TSC. Indeed, given that the voltage drop between 
left lead and TSC is V/2, linear differential conductance 
dl/dV is equal to the universal value of 2e 2 /h. Next, 
at the symmetric point = F^, Andreev current be¬ 
comes zero, and linear differential conductance between 


The frequency dependence of the cumulants for ea/T = 
—2 with different k and S are shown in Fig. 2 (left 
panel: 5 = 0; right panel: S/T = 0 . 02 ). As one can 
see, the cumulants exhibit plateaus corresponding to the 
universal values, see Eq.(49), in the frequency range 
ui < min{r, ac 2 /F} which allows one to distinguish the 
Majorana physics from the other non-Majorana effects. 
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u)/r 



(o/r 



- 0.04 - 0.02 0.00 0.02 0.04 


u)/r 


o.oo 

-0.05 



- 0.04 - 0.02 0.00 0.02 0.04 


(Jj/T 



w/r 


w/r 


FIG. 2. The cumulant spectrum C n (ui) for n = 2,3,4 for different k, 5 = 0 (left panel) and 8/T = 0.02 (right panel). Here we 
set t d /Y = —2.0, Tl = Tij. 


Next,we consider the effect of finite Majorana degen¬ 
eracy splitting 5 ^ 0 which affects the cumulant spec¬ 
trum C n (ui) at small frequencies to —> 0. One can see 
that, provided n > T, the cummulants with 5 = 0 and 
) / 0 are similar for |w| 8. Therefore, in order 

to observe the universal values of C n s, one has to ad¬ 
just the source-drain bias V to the appropriate regime: 
min{«: 2 /r,r} eV S. One can notice that there is 
also a redistribution of the spectral weight for small k/T. 
Therefore, large k/T > 1 limit is more favorable for the 
experimental measurement of the cumulants. 

Finally, we plot the second cumulant spectrum C^w) 
for ed = 0 in Fig. 3. We can see that although the quanti¬ 
tative value show small changes compared to td/T = —2 
result, the conditions for the source-drain bias shown 
above still hold indicating that our results are robust 
against changes of e d - 


IV. WEAKLY INTERACTING SPINFUL QD 
COUPLED TO A MZM 

In this section, we consider the spinful model for a QD 
coupled to a MZM which is relevant in the context of the 


Majorana proposals involving topological insulators [9, 
10, 15]. Indeed, MZM can be localized, for example, at 
the domain wall between an s-wave superconductor and 
a magnetic insulator of a Quantum Spin Hall insulator 
heterostructure. Assuming that the magnetic insulator, 
polarized along z-axis, does not affect the spin in QD (i.e. 
Zeeman splitting in QD is negligibly small), one arrives 
at the following effective Hamiltonian: 

-ffQD —mzm = ^ £ddld<j + U 

(7 

+ + dj)7i + *57172- (51) 

Here we also include the effect of inter-particle interac¬ 
tion U assuming that it is weak, i.e. U -C F, k. The 
opposite limit of strong Coulomb interaction in the dot 
is considered in Sec.V. 

As shown in the previous section, the CGF for spin¬ 
less QD becomes universal (i.e. independent of and k) 
due to the MZM coupling. In the spinful case, only one 
channel (e.g. spin-up) effectively couples to the MZM, 
see Eq.(51). Thus, the CGF will also have a nonuni- 
versal contribution from the spin-down channel which is 
decoupled from MZM. However, as follows from Eq.(35), 
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FIG. 3. The td = 0.0 result of the second cumulant spectrum C 2 {oj) for different k, 5 = 0 (left panel) and <5/r = 0.02 (right 
panel). Here we take r l = F_r. 


Normal propagator : 


Anomalous propagator : 


Coulomb interaction : 

(a) 



FIG. 4. a) Diagrammatic representation of the normal and anomalous impurity Green functions and the Coulomb interaction; 
b) and c) diagrammatic representation of the self-energy due to leading order corrections of the Coulomb interaction. 


higher-order cunrulants (i.e. n > 1) from the spin-down 
channel vanish at Cd = 0 and Tr = Tr enabling one to 
observe the universal part originating from the spin-up 
part. Thus, some fine-tuning is necessary in this case (as 
opposed to the strongly interacting case in Sec. V). In ad¬ 
dition to the aforementioned corrections to the universal 
features in FCS, one should also consider the effect of 
Coulomb interactions. Without loss of generality, we set 
Xr = 0 and calculate effect interactions on charge fluctu¬ 
ations through the left lead. Our conclusions also apply 
to charge fluctuations through the right lead. 

We now consider effect of weak interactions U -C 
{r, k} on FCS. We first calculate the contribution of 
Coulomb interaction to the self-energy and then obtain 
the corrections to the CGF in powers of U. Up to the 
second order in U, the corresponding Feynman diagrams 
are shown in Fig. 4. The linear in U contribution to 
the self-energy E A |, see Fig. 4, merely represents the 
renormalization of the QD energy level ed- This is a triv¬ 
ial interaction effect which does not modify our previous 
conclusions. We, therefore, focus on U 2 contributions. 
The self-energy in the Narnbu space has the following 
form 


S At = 




(52) 


where E^'' and E^J are particle-hole conjugation of E A | 

and E^d- We note that all the Green’s functions here 
depend on the counting field \r. The details of the cal¬ 
culation of E Ai is presented in the Appendix B 1. After 
some manipulations, the cumulant generating function 
for each spin-channel can be written as 




det 


[Qdd,U =c 


i -i 


- E At 


det 


rn A i=° 1 1 _ vAi.=o 

l^dd,U=0i ^ 


1 r 

2 /_„ 


dw, det 

\Qdd,U= o] 

27r n det 

[ o X l=0 1 _1 

L^dd,!7=0j 


(53) 


r 


dui 


det 


In 


^4x4 Qdlu^ 


-oo 27T (Jgt 


14X4 ^ Qdd(u=0 EAi=0 


where 


[Qdd,u= o] _1 = [Qtid] _1 - E ( E M T, a ® iV&.a 

a. i 

x(E M Ua®y)> (54) 


The first term in the generating function corresponds to 
the result for noninteracting case, see Eq. (30) whereas 
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the second term originates from the interaction-induced 
corrections. We consider weak Coulomb interactions 
U {T, A, eF}, and keep the leading order terms in U: 
Qdd,u=o'^- ,XL ~ C^ 2 /min{r, A, eV } 2 < 1. By expanding 
the second term in Eq.(53) up to quadratic order in U, 
we obtain 


lnx(A) 


1- 

2 ^ 27t 


det 

7 

o' 

II 

b 

<< "8 

det 

a- ^ 
a-h 

II ° 
o 

1 


r 


dco r 
27T - 


Tr (Qdd,u=o sAL ) 


— Tr (Qdd ,u=o ^ Al = °) 


(55) 


where we used the relation det (I + xA) ~ 1 + xTr(A) for 

x -C 1. 

After some manipulations (see Appendix B 1), the cu- 
mulant generating function for small U can be written 
as 


ln.Xix(Ap, U ) 


In \cr(^L,U = 0) 



(56) 


Note that the matrix formalism of the Green functions 
has the form ( ) (he. with a L-0 rotation). How¬ 

ever, for the convenience of the calculation, we will 
also use the matrix Green functions in the Schwinger - 
Keldysh space (without L-0 rotation). Here the function 
E\ (A here means (Ap,Ap)) in the Schwinger - Keldysh 
space has the form 


^A — 






(n? lt (-fi)n^(fi) + < t (-H)n^(H)), 
^ + n| T (-0)n^(f))), 

fj (n< t (-H)H< 4 (H) + H< t (-H)H< 4 (H)), 
g (n> T (-H)H> 4 (H) + H> f (-0)H> 4 (H)), 

(58) 


and 


CT if 

“A 



(59) 


The polarization functions Ilp^ and Ilp CT (see Fig. 4) 
are calculated in Appendix Bl. One can notice that, 
for ed = 0, the particle parts are exactly the same as 
the hole parts: np i(T = np jCr . Thus, the functions 
for spin-up channel have the same form as those for spin- 
down channel, and we drop the spin-index in Ed from now 
on. Lastly, it is well-known that the function = 0 

(here A = 0 means Ap = Ap = 0) vanishes because of 
causality and unitarity. However, the presence of artifi¬ 
cial counting field A (t) having different sign on forward 
branch and backward branches of the Keldysh contour 
breaks unitarity. Therefore, the relation S^l 0 = 0 does 
not hold for A A 0, and we have to evaluate it explicitly. 
In general, the calculation of 5^"^ is not very illuminat¬ 
ing but some simplification can be obtained by expand¬ 
ing the interaction-induced corrections in powers of eV 
assuming that eV 0. After some manupilations (see 
Appendices B2, B3, and B4), we find that the leading 
contribution to CGF is proportional to V 3 : 


—if 
^A ,A 


(eV) 3 ( 


„A\(3) 
A ,A 


“A ,B J ) 


(60) 


where the functions S 


if,(3) 
A, A 


and 5 


if,(3) 
A ,B 


are given by 


where 


.if,( 3 ) _ {e~ iXh - 1) ^ , , o tt% (e iX - - 1) [(e iXL - 1) - 2(e iA - + 3)« 2 ] 
j a,al — 6?r p4 a ) 


?K,( 3) 
J A .B 


1 


247r 3 T 4 K 2 (l + e iXL Y 
+ (T 2 - 8 K 2 )(e~ 2iXL - 1) 


4 ' 127r 3 (e lAi + 1) 2 T 4 k 2 

15n 2 (e iXL - 1) + K 2 (e 2iXL - l) + (—2T 2 - 15K 2 )(e- iAi - l) 


np2, t ( 0+ .°)-np2, t (o + ,o) / _ 


P 2,tV 


67T 2 r 3 


[e~ iXL — l), 


(61) 


(62) 


O («) = 


dfl 

27T 


n£ t (-n,o)^=-+ n£ t (-n,o)^=- 

(|H| + i) 2 p ’ n (\n\ ^*) 2 


~ „ K 

, with S l = K = —. 

. r > r 


(63) 


One can check that indeed above expression vanish for Ap = 0, as required by unitarity. The integral in Eq. (61) is a 
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FIG. 5. Diagrammatic representation of the cumulant generating function up to the leading order corrections of the Coulomb 
interaction. 


dimensionless number depending only on n/T. The functions np 2 T ^(0 + , 0) in Eqs. (62) are defined in Eq.(B19); these 
functions depends on T and k. [88-91] [92] 

Combining all the terms, the final expression for the cumulant generating function of the spinful QD at 0{U 2 )-level 
can be written as 


lnx(A L ,£7) « ]n X t(X L ,U = 0 ) + ln^(A L ,U = 0) - TU\eVf(^lf 1 + sj£ 3) ). 
and the leading correction for the current and shot noise (for left junction current) can be written as 

u = e u 2 {eVf f Q(Ti) _ 1 ng 2 , t (0 + , 0) - n£ 2 , t (0+ 0) \ 

L h T 4 Y 67T \ 4 J 6-7T 3 27T 3 67T 2 J 1 

u e U 2 (eV) 3 f Q(«) / ^ 2 \ 1 + 2 k 2 1 - 62 k 2 n £ 2 T ( 0 +, 0 ) ~ n £ 2 t ( 0 +, 0 ) 

LL h E 4 y 67r 4 J 247t 3 k 2 487t 3 k 2 67t 2 


(64) 


(65) 

( 66 ) 


Thus, the leading order correction to the generating function in the presence of MZM coupling is of the order of 
(eV) 3 which is the same in the case without MZM considered in Ref. 68. Since the leading order correction to the 
cumulants C n is of the order of {eV) 2 , the interaction-induced corrections do not affect cumulants at small bias. 
Therefore, we expect that one can observe the universal values of the cumulants discussed in Sec. Ill in realistic 
experimental conditions. 


V. STRONGLY INTERACTING SPINFUL QD 
WITH MZM: INTERPLAY OF KONDO AND 
MAJORANA COUPLINGS 

A. Theoretical model 

In this section, we study another nontrivial case cor¬ 
responding to a strongly interacting spinful QD coupled 
to a MZM. The Hamiltonian for the system reads 

TIqd-mzm = E Cddlda + Ud^d^d^di 

G 

+in(d^ + dj ) 7 i + *< 57172 - ( 67 ) 

Once again, we assume here that Zeeman splitting in a 
QD is negligibly small (see discussion after Eq.(51)). For 
large on-site Coulomb interaction U and single-electron 
occupancy, we have to consider interplay of Kondo and 
Majorana physics [46]. In the limit of single-electron oc¬ 
cupancy {T, k} <C |e<j| < U, one can study the prob¬ 
lem using a slave boson mean field (SBMF) approxi¬ 
mation originally developed for an infinite-C/ Anderson 
model [93, 94]. This approach allows one to eliminate 
double occupancy in the QD and significantly simplify 
the problem. 


For the sake of completeness, we outline here main 
steps of SBMF approach. We refer a reader to Ref. 46 
for more details of SBMF calculation in the presence of 
MZM. We first rewrite fermion operators in QD in terms 
of the auxiliary boson b and fermion f a operators, i.e. 
da- —t fab - This procedure requires to introduce a con¬ 
straint tfb + J2a faf° = 1 on the Hilbert space. After 
the transformation, effective Hamiltonian becomes 


#SBMF = i^Leads + E £ dfif° + *«7i(/t^ + fib) 

G 

+ E E tai'f’i^Wfab^ + h.c.) + 7(57172 ( 68 ) 

CX—L.R G 


where the lead Hamiltonian I?Leads is unchanged. Next, 
we apply mean field approximation and replace the 
bosonic operator b and the Lagrangian multiplier 77 en¬ 
forcing the constraint by their mean-field expectation val¬ 
ues. We choose ( b) = (b') = b to be a real positive num¬ 
ber. The mean field parameter b and 77 can be determined 
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Kir 





FIG. 6. Dependence of the cumulants C n ( 0) = {{S n q))/M (5 = 0 and eV —> 0 limit) on the Majorana coupling k in Kondo 
regime from SBMF calculation. We choose Yl td/Y = —10.0, and band width A/F = 30.0. 



FIG. 7. The cumulants C n (0) for n = 1,2,3,4 as a function k and 5. Here td/Y = —10.0, Yl = r_R, eV/Y = 0.001, and 
A/r = 30.0. 


self-consistently by minimizing the free energy [46]: 

b 2 + J2(fUa) = l, (69) 

a 

26?y +1 E E((/^,a(0))+C.C.) 

a—L,R a 

+in(li (fl + ft)) = 0 - 


Here, we assume the eV max \ T K . k} and thus neglect 
the dependence on voltage bias eV in the SMBF calcu¬ 
lations, see also discussion in Ref. [95]. Above equations 
determine thermodynamics of the system, and we are 
now ready to compute transport properties. Once the 
mean field values of the auxiliary parameters are deter¬ 
mined in SBMF approximation, spin-up and spin-down 


(70) 
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channels become decoupled. Thus, effectively the prob¬ 
lem reduces to the previous model - spin-up channel is 
coupled to a MZM whereas the spin-down channel is not. 
The chemical potential and couplings are renormalized: 

~t + ih t a —> t a = bt a , r Q —> T a = b 2 T a , and 

k —>■ k = bn. 

Using the results of Sec. Ill, it is rather straightforward 
to obtain the CGF at small voltage bias eV —> 0 in this 
case 


In ,\km 
M 


= In 

eV —>-0 


/ r L e 1 ^ + r R e-* x « \ 

V r L + f R J 


+ In 



4r L r^ 

+ (f L + f *) 2 


^gi(Ai-AR) 


1) (71) 


Here the first and second terms correspond to the spin- 
up and spin-down channels, respectively. As discussed 
below, the renormalized QD energy is close to the 
Fermi level, i.e. ej, —> 0, for a large parameter range 
{/■c,r} < \ed\ (so-called universal limit). Thus, in the 
case of symmetric right-left lead couplings = F i; the 
second term in Eq. (71) is simply given by i(X l — Xr) 
and, therefore, does not contribute beyond the first cu- 
mulant. Thus, the shot noise as well as other higher order 
cumulants are once again given by the universal pattern: 


{Ci(0) • • • C 8 (0)} = 


3 1 
2’ 4 


, 0 ,- 4 , 0 , 4 .°,— 


17 

16 


C n ( 0) = for n > 1. 


(72) 


with E n {x ) being the Euler polynomial. 


B. Results and Discussion 

We now analyze different cases in details. We be¬ 
gin with the case of zero degeneracy splitting <5 = 0. 
The dependence of the cumulants on n for <5 = 0 and 
eV —>■ 0 is numerically shown in Fig. 6. A recent 
study based on SBMF approach [46] shows that there 
is a crossover between Kondo- and Majorana-dominated 
regimes as a functi on of the MZM coupling k. For 
k < k c = y/T K /T\e d \, the mean field solution is de¬ 
termined by the Kondo temperature Tr' 

f = F6 2 =T x = Aexp(-7r|e d |/2r) (73) 

and the renormalized energy level is id = \td + ry| ~ 
Tb 4 . Here A is the bandwidth and T = F# + .JSince 
b <C 1, the renormalized energy ed is small <C F [46]. 
Thus, the cumulants are given by Eq.(72). In the case 
of intermediate MZM coupling k k c , the parameter 
b ~ n/\ed\ is determined by the Majorana coupling rather 
than the Kondo temperature. Still, however, if |ed| ^ 
r, re, the occupation of the empty state is small b -C 1, 
and the position of the renormalized level is still close to 
the Fermi energy Q ~ Tb 4 <C T [46]. Thus, the cumulants 
are given by Eq.(72) for T# = F^. 


Next we consider the effect of a finite energy splitting 
6 ^ 0 and voltage bias on our prediction, which is impor¬ 
tant for the experimental detection in realistic settings. 
The cumulants C„(0) for n = 1,2,3,4 as a function k 
and S are shown in Fig. 7, where we focus on the limit 
\ed\ {F, k}. One can see that in order to resolve the 
universal quantized values, one has to adjust the voltage 
bias in the following range min{K 2 /r, F6 2 } eV <5 
where b = ^T K /T for n <C k c and b = «/1 | for 

k c •C k •C |ed|- The plot of the cumulant power spectra 
as a function of the splitting 5 and Majorana coupling k 
is shown in Fig. 7. One can notice that the width of the 
plateau around the quantized values gradually shrinks 
with increasing 6. 

Finally, we consider strong MZM coupling limit k > 
\e d \ such that b ~ 1 and ~ T. In this case, although 
the energy level shift does not affect the universal values 
for spin-up channel (due to MZM coupling), the posi¬ 
tion of the renormalized level q affects the cumulants in 
the spin-down channel. Therefore, overall results become 
nonuniversal and depend on microscopic details such ed 
and k. As Majorana coupling k is increasing, cumulants 
start to deviate from the universal values as shown in 
Fig. 6. 


VI. CONCLUSIONS 

In this paper, we study the full counting statistics of 
charge fluctuations in a QD device with a side-coupled 
TSC as shown in Fig.l. Two normal metal leads, cou¬ 
pled to the QD, are also introduced in order to detect the 
charge tunneling events. Using a Keldysh path-integral 
approach, we compute the cumulant generating function 
for the QD with MZM coupling in this two lead struc¬ 
ture. We first consider a noninteracting spinless system 
and find that for the symmetric left-right lead couplings, 
the zero-frequency cumulants exhibit a universal pattern 
described by a series of numbers generated from the Euler 
polynomial. This result is independent of the microscopic 
parameters of the dot, i.e. QD energy level and QD-MZM 
coupling. For small topological degeneracy splitting due 
to a finite-size TSC, the universal pattern can also be ap¬ 
proximately observed provided the voltage bias is in the 
appropriate range. We also compute FCS for a spinful 
QD setup, and consider effect of Coulomb interactions 
both in the perturbative (weak interactions) and non- 
perturbative (strong interactions) regimes. In the former 
case, we find that the interaction-induced corrections to 
the cumulants appear in the cubic order in voltage V and 
quadratic order in interactions U indicating that the uni¬ 
versal pattern characteristic to Majorana zero modes can 
be, in principle, measured experimentally. We find that 
the optimal regime for that corresponds to the strongly- 
interacting spinful QD in the single-occupancy regime 
with the intermediate MZM coupling. Our results pro¬ 
vide a complete tool set for detecting Majorana modes in 
tunneling transport measurements which goes far beyond 
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Appendix A: Derivation of the coefficients in CGF Eq. (37) 


In this Appendix, we provide details of the calculations of FCS for non-interacting case. To evaluate the generating 
function, we first define 


K(X L ,X R ) = det 14 x 4 - Qo : dd^2 (^2 M T,a ® 7*) t( 3o,a(^ M T,a ® 7*) 


(Al) 


We consider a symmetric source-drain bias {hl = eV/2 and n R = —eV/2), and insert Eq. (19), (22), and (27), (28) 
into K(Ai,Aij), and obtain 

K(X L = o, Afl = 0) = 1 + (T L + r R ) 2 ([G R dd ] 2 + (r L + r R ) 2 [F R dd }* - 2 [F 0 * rf ] 2 ( -1 + (r L + r R ) 2 G R dd G R dd ) 

+K J 2 (l + (I# + T R ) 2 [G R dd ] 2 )) (A2) 


and 


K(Al, X r ) = K(0, 0) - Citil(1 - n L ) - C 2 7i_r(1 - n R ) + Binz,(l - n R ) + B 2 71r(1 - n L ) 

+Wn L n R (l - n L )[l - nR) + In L n R (n L - n R ), (A3) 


where 


Cl = 4e~ i(Al ' +AR) r z ,r ii ^4T L r“ A h) (e iXL - e iX «) 2 ([F R dd ] 2 - G R dd G R dd f - (1 - e^ +x ^) 2 [F R dd ] 2 ) , (A4) 
c 2 = 4 e - i ( Ai + A «)r L r ii ^4r L r fl e _i(Ai_AR) (e iAi - # Ar ) 2 ([E 0 * d ] 2 - G R dd G R dd f - (1 - e^ +x ^) 2 [F 0 R dd } 2 ) , (A5) 

Bi = 4e -2iAi (-2e’ AR (# Ar - # Ar ) ([F R dd ] 2 - G R dd G R dd ) 2 T 3 L T R + e 2iA - ( - 1 + e 2lA «) [F R dd ] 2 T 2 R 
^ {e iXh - # Ar ) T l T r ([G« J 2 + [G R dd ] 2 + 2 ([F R dd ] 2 G R dd G R dd ) " T 2 ^ 

+ r l (- (-1 + ? 2iXl ) K dd ] 2 + 4 # Ar (—e iAi + # Ar ) {[F R dd ] 2 - G R dd G R dd )" F^ ), (A6) 

B 2 = 4e -2iAi ( - 2# Ar (e iAR - # Ar ) {[F R dd ] 2 G R dd G R dd )" r£r* 

+e 2lAi (-1 + e 2iAR ) [F R dd ] 2 T 2 R (# Ar - # Ar ) r L r* ([G R dd ] 2 + [G R dd ] 2 + 2 (\F R dd } 2 G R dd G R dd ) " T 2 ^ 

+ r l (- (-1 + e 2iXL ) [. F R dd ] 2 + 4e lAR (-e iA - + # Ar ) ([F* dd } 2 - G R dd G R dd )" T 2 ^ ), (A7) 

F = 256r 2 r^([F£j 2 - G R dd G R dd ) 2 (sin A ^ ii ) 4 , (A8) 

1 = 128*r 2 T^ {[F R dd ] 2 - G R dd G R dd ) 2 ( sin * sin(A L - X R ). (A9) 

where the Green functions, i.e. G o dd ^ Gq ,ddi and Fo jdd , are defined in Eq. (23) (24) (25) of the main text. 
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Appendix B: Calculation of the CGF corrections due to weak interaction effects 


1. Derivation of the interaction correction formula 


For the convenience of the calculation, we consider the polarization function ( see Fig. 4 ) in the Schwinger - 
Keldysh space (without Larkin-Ovchinnikov rotation) 


ftp (ft) 


fi<(n)\ 
\n>(n) n£(fi)y 


r°° duJi_ + f2 ) G 'L“( w i) + F Jd(ui +tyFjj{ui) 

J-oo 2vr l v G>> 1 + n)G<> 1 ) + F>(u; 1 +n)F^(a; 1 ) 


G dd( ui +i)+ F dd( w i+/ Bn 

G I>i + + i&(wi + n)i^( W i) y j 


Here the subscript P indicates the particle channel, its particle-hole conjugation (i.e. the hole channel) IIh(H) has 
the same form but with replacement G dd —► G dd and F dd —> F dd ■ The respective self-energy for the spin-up channel 
can be extracted from 


^, t M 


iU 2 


iU 2 


iU 2 


iU 2 


dn 

dn 

2?r 

dn (Fj d ^-n) n£ 4 ( 0 ) 

dnfc^-n^n) 

27r 


G ds^-mi^)\ 

Fj St (co-n) nUp(fi)J 

Fdd.,^ ~ ^)nH,4.(^)\ 

G L, t 


(B2) 

(B3) 

(B4) 

(B5) 


where the spin-up channel couples to MZM and the spin-down channel does not. We want to calculate the following 
function, i.e. interaction correction in Eq.(55) 


/l£ Tr ( Qdd, t ,u=oH^ L M) 


duj 

2tt 

dui 


Tr(G da;t ME^ t H + F ddit (u;)E^ )t ( w ) + F dd|t ( W )E^ f (u,) + G ddit ME^ f (u,) 


^Tr( 7 cl G ddiT ( W ) 7 cl E^ jt ( w) + 7 cl F ddjT (w) 7 cl E)T jt (w) + 7 cl F ddiT M 7 cl S^, T M + 7 cl G ddit M 7 cl E^ t M 

(B6) 

where 7 d = 12x2- Note that the matrix Green functions and self-energies (which are from Eq.(55)) have the form 
( qk ( qa ) (he. with the L-0 rotation). Therefore, the function above just corresponds to the classical-classical part 
( K ) of a certain polarization function, and we define a function for interaction correction 


TT 2 'ZlK 

u —A,t 


^Tr( 7 cl G dd>t (o;) 7 cl E^ it ( W )+ 7 cl F dd7 (u;) 7 cl E^ it ( W )+ 7 cl F ddit (u;) 7 cl E^ it (u;)+ 7 cl G ddjt ( W ) 7 cl S^ iT ( W 

. ( B7 ) 

which can be described by the diagrams shown in Fig. 5. The function S;, in the Schwinger - Keldysh space (without 
L-0 rotation) has the simple form 






(B8) 
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where 


and 



(n?, t (-n)n£ 4 (fi) + n£ it (-n)n£ 4 (n)), 
(nf it (-n)ft? 4 (n) + n^ iT (-n)nj 4 (n)), 
g (n< t (-n)n< 4 (fi) + n< t (-n)n< 4 (n)), 
(n£, t (-n)n> 4 (n) + n> t (-n)n> 4 (n)), 

sf=(sr + sf- 3 <-3>)/2. 


Finally, we reach the formula in Eq.(56) of the main text. 


(B9) 

(BIO) 

(BH) 

(B12) 

(B13) 


2. Expansion of 

Now, let’s look at how the correction changes as 
a function of source-drain bias eV for T = 0. For small 
eV, we can expand the function 

Sf (eV) = Hf (0) + (eE)sf’ (1) (0) + (eH) 2 ~f' (2) (0) 
+(eH) 3 sf’ (3) (0) + --- (B14) 

Due to the causality reasons (this is still true in the pres¬ 
ence of artificial counting held), the polarization func¬ 
tion at T = 0 and eV = 0 has the following properties: 
np CT (f2) oc n( fi) and rip CT (fi) cx 1 — n(Q), where n(fl ) = 
1 — 9(fl) is the Fermi-distribution function at T = 0. 
Then, we find Ef(eV = 0) = 0 and E^(eV = 0) = 0 
for T = 0. Although the time-ordered and anti-time- 

ordered parts of the Green functions (Gjj T ,G^ T , F^’ T , 

and ) depends on the counting held Ap, their de¬ 
pendence on Ap enters in a way such that all the A l 
dependent terms have a pre-factor np( 1 — nit). At T = 0 
and eV = 0, np( 1 — n R ) vanishes and time-order and 
anti-time-order Green functions are independent of Ap. 
Therefore, we have E^(eV = 0) = S^ =0 (eV = 0) and 
S J)(eV = 0) = E^ =0 (eV = 0). Combining those rela¬ 
tions with Sf =0 (eV = 0) = Ej) =0 (eV = 0) + Sj =0 (eF = 
0) = 0, we can show 


sf (0) = 0. (B15) 

This is intuitively obvious since there should be no cur¬ 
rent in equilibrium (i.e. eV —> 0). 

To consider higher order terms, we expand the inte¬ 
grand of the polarization function integral in order of 
eV. These integrands have the following form 

F a (u)i,n\nL (wi + £l),n L (ui),n R (u!i + tt),n R (u>i)) 

= G^( Ul + D)G“ d -( Wl ) + F£(wi + ^fd(wi)(B 16 ) 


where a = T,T,<,>, The bias eV only enters through 
the Fermi distribution function n k and n R : ul(uj) = 
1 — 0(w — eV/2), rifl(w) = 1 — Q(u; + eV/2). Due to the 
properties of Heaviside theta function, by expansion and 
resummation, one can prove the following relation: If 
one has a series of functions ni(u>) = 1 — 9(uj — uq) with 
* = 1,2,--- ,k and uq ^ lo 2 ^ ^ then 


n 1 


F(m,n 2 ,--- ,n k ) = 

F( 0,0,-- - ,0) + [f(1,1,--- , 1) — F(0, 1, ■ ■ ■ , 1) 

i i +1 

F(oT^O,1,--- , 1) - F(o7^To, 1, ■ • • ,1)U 


F{0, ■ ■ ■ , 0,1) — F(0,0, • • • , 0) 


n k . 


(B17) 


Note that this formula depends on the order of the ar¬ 
guments in the Fermi distribution function. We then de¬ 
fine the polarization function for different regions where 
the function F has different forms. First of all, we con¬ 
sider > eV such that —eV/2 — D < eV/2 — Q < 
—eV/2 < eV/2, and define 


fl a P 1 jn,eV)=i 



doji 

2n 


(wi,n), 


(BIB) 


and Ff (uq, fi | n R (uq + fi), n L (uq + H), n R (uq ), n L (uq ) ) = 
F a (uq, H| n L (uq + Q.) , n L (uq ), n R (uq + D), n R (uq )), where 
F a is the same function as the one defined in Eq. (B16), 
but the arguments in the function F“ have the different 
order. Secondly, we consider 0 < H < eV such that 
—eH/2 — f2 < —eV/2 < eV/2 — D < eV/2, and define 

/ OO 7 

(B19) 

where F.f(uq, fi| n R (uq + f2), n R {ui\), 7ip(uq + 
fi),nz,(uq)) = F a (uq,D|ni(uq + D),n L (uq),n fl (uq + 
H),n R (oo 1 )). Thirdly, we consider —eV < Q < 0 such 
that -eV/2 < -eV/2 -Q < eV/2 < eV/2 - D, and 
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define 

/ OO 7 

(B20) 

where -F 3 a (wi,fl|n fl (a;i),n^Wi + fi),« l(wi), n p (wi + 
fl)) = F a (wi,fi|nL(wi+fi),n p (wi),n p (u;i + fi),n p (wi)). 
Finally, we consider fi < —eV such that —eV/2 < 


Following the definition above, the correction ^ can bi 

^o; _ ^a. 

~ “A, 


eV/2 < —eV/2 — fi < eV/2 — fi, and define 

/ OO 7 

^F“K,fi), (B21) 

where -F“(wi, fi|n p (u;i), n p (wi), n p (wi + fi),n p (wi + 
fi)) = F a (wi,fi|nL(wi + fi),n i (wi),n fi (o;i + fi),n_ R (wi)). 

written as (here we consider e,j = 0) 

+ S“ B , (B22) 


where 


■"A.A 


J A ,B 


'°° dfi - . f° dfi - 

= 2 / —n? lit (fi)n P 44 (-fi)“ + 2 / —n£ 4 it (fi)n% u (-fi) 


'0 


r eV 


2n 
dSl r 


' —OO 


n£ 2 jt (fi)n P 3 ii (-fi)“ - n^ 1 )t (fi)n P 4 , i (-fi) c 


-Vo 2 tt L 

/■° r/fi r - 

+2 / — n^ 3 it (fi)n P 2 i 4.(-fi) Q - n^ 4 ;t (fi)n Plji (-fi) c 


' —eV 


2tt . 


r-’CK 

'A.A 


'-70: 

-A,S- 


(B23) 


Here the factor 2 comes from the summation of both particle and hole part (note that li Po . = n P , CT for = 0). 
The whole correction includes two parts: the first part and a leftover part S“ B . To check this formalism, we 
considered the case without MZM coupling and reproduce the result shown in Eq. (21) of Gogolin and Komnik [68]. 
Note, in order to obtain the right result, we need to take appropriate order of limits. If we want to reproduce the 
k = 0 result, we have to take the limit k —>■ 0 before taking the w —» 0 limit for linear response. Then, we focus on 
the case with MZM coupling (where if we want to consider zero Majorana splitting, i.e. <5 = 0, we have to take the 
limit S —> 0 first before taking the ui —> 0 limit for linear response). 


3. Derivation of the first part of interaction correction E*^ 2 ' 


First of all, we focus on the first part A , and expand it in order of eV.We also notice that after the transformation 
of Eq. (B17), the integrands of the polarization functions are linear in the Fermi distribution function. Therefore, 
the expansion can be obtained analytically by expanding the Fermi distribution function 

eV / eV\ 1 ’ / eV\ 2 1 " / eV\ 3 

n{u±— )=n( W )-%)(± T ) --dM(±—) +••• (B24) 

After expanding the Fermi distribution function in order of eV, we further expand the polarization function 


n P (fi, eV) = n P (fi, 0) + (eE)flp ) (fi, 0) + (eH) 2 ri^ 2) (fi, 0) + (eE) 3 n^ 3) (fi, 0) + • • • 


After the integration by part for the Dirac-delta function, the linear terms can be obtained 


n“’ (1) (fi > 0,0) 

llp’ (1) (fi < 0,0) 


--iT(-fi, fill, 1,1,1) + 2F“(—fi, fi|0,1,1,1) - Ff(- fi, fi|0,0,1,1) 

47r L 

-iT( 0, fi|0,0,1,1) + 2E“(0, fi|0,0,0,1) - Ff(0, fi|0,0,0,0)], 

V [ - F?( 0, fi|l, 1, 1, 1) + 2F“(0, fi|0,1, 1, 1) - E?(0, fi|0,0,1,1) 
-F7(-fi, fi|0,0,1,1) + 2F 4 “(-fi, fi|0,0,0,1) - F“(—fi, fi|0,0,0,0)1, 


(B25) 


(B26) 

(B27) 


where a = T, T, <, >. After simplification of those expression above, we find for both spin-channels at T = 0 

rip (1) (fi,0) =0, flp (1) (fi,0) = 0, flp (1) (fi,0) ~ 0(-fi), rip’ (1) (fi,0) ~ 0(fi). (B28) 
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We then obtain the linear correction to the generating function 

sf (1) (o) = 2 f ^[n^-fi^n^fi^ 

+n< t (—Ci, o)n^ 1} (ci, o) + n> t (-ci, o)ii>f\ci, o) 

+ftp^ 1) (-n, o)ft£ 4 (fi, o) + o)n? 4 (fi, o) 

+tipf- ] (-ci, o)n< 4 (fi, o) + } (-fi, o)n> 4 (fi, o) 

= 0. (B29) 

Similarly, the quadratic terms can be written as 

ri“’ (2) (n > o,o) = = -n,n|i, i,i,i) -d ul Ff( Ul = -fi,fi|o,o,i,i) 

1 D 7 T 

+a U)1 F 1 a (w 1 = o, fi|o, o, i, i) - d ul F?(u! = o, fi|o, o, o, o)], (B30) 

n“’ (2) (n < o,o) = = -n,n|i, i,i,i) -d ul F^( Ul = -fi,fi|o,o,i,i) 

1 D 7 T 

+a wi F“(wi = 0,O|0,0,1,1) -d Ul F^( Ul = 0,0|0,0,0,0)]. (B31) 

The quadratic term of the generating function is therefore 

sf (2) (o) = 2 f ^[n^-fi^n^f^^ 

+Upf\-n,o)ul l (n 1 o) + np’{ 2) (-fi,o)nf; 4 (fi,o)l, (B32) 

where we use the relations fl > (— Cl, 0)fl > (n, 0) = 0, II < (— Cl, 0)fl < (Cl, 0) = 0, n^ 1 ^’ T (—fl, 0) = 0, and n(TT(_ Cl, 0) = 
0. After the simplification of polarization functions , we find that the time-ordered and anti-time ordered parts 
f[p’^(fi,0) (for a = T,T) do not depend on the counting field A l- Then, we conclude that the quadratic term 
^.a.U)( 0 ) j g game ag term w jth Al = 0, which is zero due to causality and unitarity 

sf (2) (0)=0. (B33) 

The cubic term of the polarization function corresponds to the integral of 5 , and thus reads 


n“’ (3) (fl > 0,0) = 


- dl F?(-Cl, fi|l, 1,1,1) + 2dlF?(-Cl, fi|0,1,1,1) - dlF?(-Cl, fi|0,0,1,1) 


^ ^ ^ 3!2327 r |_ Wuj i ^ i ^ > i > > ’ > ' u>i i- v “i 

-^^“(O, «|0,0,1,1) + 25^(0, fi|0,0,1,1) - dl x Ff( 0, fi|0,0,0,0) 


np (3) (fi < 0,0) = [ - dlF%( 0, fill, 1,1,1) + 2dlF2(0, fi|0,1,1,1) - 3^(0, fi|0,0,1,1) 

-d 2 Ul F2(-Cl, fi|0,0,1,1) + 2dl 1 F2(-Cl, fi|0,0,1,1) - ^Ffi-Cl, fi|0,0,0,0)]. (B35) 

We evaluate and simplify the function above for td = 0 and Tl = Tp = Y/2, and obtain 
-'Tis) , —i(e iXL — l)\(e iXL — 1)T 2 — 2(e iXh + 3)k 2 1 1 

n 5t } ( n .°) = —- /'..a. , o aT J 2 - — no. ■ ( B36 ) 


flpf\ci,0) = 
np;f(fi,o) = 
nJ{ 3) (fi,o) = 


Tp = T/2, and obtain 

2(e iXL + 3 )k 2 ] 

1 

2 k 2 

(|fi|+7T) 2 

2(e iXL + 3 )k 2 ] 

1 

2 k 2 

(|fi| - iT) 2 


—i(e iXl — 1) 1 

12t rT 2 (|fi| +iT) 2 ’ 

—i(e~ lXL — 1) 1 

127tT 2 (|fi| -iT) 2 ' 
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The cubic term of the generating function reads 


~ A ’ (3) (0) = 2 


U X,A 


r°° dn 

— on 27T - 


n£ it (-n, o)fip;} 3) (fi, o) + o)ft£{ 3) (fi, o) 

T T >(3)/_n mffT fo ro tt'^’T 3 )/ 


+np;| JJ (-o, o)n^(n, o) + n P y’(-n, o)n^(n, o) 
- <e - A ‘- 1) rg[ft?, t( -n,o ); 1 


67TT 2 7-00 2ttL a ’ tv ’ ; (|fl|+«r) 2 

- i(e^ - 1) [(e iAij - l)r 2 - 2(e* Ai + 3)« 2 ] 


+ ri£ t (-n,o) 


i 


m-iTY 


247r(e iAi + 1) 2 T 2 «; 2 


(e 


— oo 
—iAi 


27T - 
- 1 ) 


n?,i(- n .o) 


i 


+ npi(— n, o)- 


(|fi| + *r) 2 1 ““’~'(|fi|-ir) 2 J 


tt 2 (e^ - 1) [(e iAi - 1) - 2(e iAi + 3 )k 2 } 
GttT 4 W + l 4 127r 3 (e iAi + 1) 2 T 4 k 2 


(B40) 


where 


o(«)= / g[n? tf (-n,o) 


(M + *) 2 


+ n£ t (-n,o) 


(M -*) 2 


, with n = —, k = —. (B41) 


This correction is nonzero for Ap ^ 0. 


4. Derivation of the second part of interaction correction 


J A ,B 


Now, let’s consider the leftover terms 


r eV 


/ o 


d,n ' 

2n . 


np 2it (fi,eV)n^(-fi,eV) - n£ ljt (n, eV)n? 4)4 (-n, eV) 


r° rfOr. 

+2 / — n^ 3T (f],eV r )n?, 24 (-f 2 ,eV r ) -n% 4iT (n,ey)n^ u (-n,ey) 


/- e y 

/■ e y 


27T . 

dn' 


4 ./n 2-7T L 


np 2 ,t( n . eV)n% 34 (-n, eV) - n£ lit (n, eV)n£ 44 (-n, e v) 


(B42) 


We want to expand the following integral in order of eV 


reV 


dn 


'o 


—n a PiA (n,eV)u a PU (-n,eV) 


i 

27T 1 
1 

27T ' 


= ^( e y)n^ >t (o+,ey)n^. 4 (o-,ey) 

„ 1 ,dtl%J0+,eV) „ , ^ an^.fO-.eV) 

■ {eV) y.(-dn -- n mr.^)- n p I ,t(0 + .^)^— z 

1 dii“ Pi!t (n,eV)ii« Pjti (-n,eV) 


+ 2^3! 


on 


O-S-0+ 


(B43) 


In the next step, we will expand the functions lIp i)(T ( 0 ± ,eV) dnflp ia (0 ± ,eV) and <9prip ijCT (0 ± ,eV’) in the order of 
eV. The way to expand ftp, a (O 1 * 1 , eV) can be found in appendix B 3. For the daflp i eV), for example, we notice 
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that 

dfl a P2 (^eV) 

on 


i 8 r 

— — / F 2 a (wi,Sl|n,R(wi +fl),n R (u 1 ),n L (u}i + n),n L (ui))dui 
^ on J- oo 


s 


{d n F 2 (io 1: fi|0,0,0,0) + [SnF 2 ( Wl , 0|1,1,1,1) - d n F 2 {u u n|0,1, 1,1 )]n fl ( Wl + fi) 


+ [ftii ? 2 (wi ) n|0, 1,1,1) - d n F 2 {L 0 1 , fi|0,0,1,1)] n*(wi) 

+ [dnFifa, O|0,0,1,1) - dnF 2 {cu u fi|0,0,0,1)] n L (u i + fi) 

+ [d n F 2 ( Wl , fi|0,0,0,1) - dnF 2 {co i, fi|0,0,0,0)]?z L ( W i)}dwi 


F 2 (-n-—,n\i, 1 , 1 , 1 ) 

eV 

F 2 (—Q + —, f2|0,0,1,1) 


,Q|0,1,1,1) 

eV 

F 2 (—fl + —, H|0,0,01,1) 


(B44) 


Similarly, for eV) with j = 1, 3,4. Then, to expand the whole function, we just need to expand the Fermi 

distribution function and the function F 2 in order of eV, which is very straightforward. For the second derivative 
CT (0 , eV), we can apply the same strategy. Combing all the terms together, we recover Eq. (62). 
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